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GENERALIZED LINEAR COVARIANCE ANALYSIS 

J. Russell Carpenter and F. Landis Markleyt 


We review and extend in two directions the results of prior work on generalized covariance 
analysis methods. This prior work allowed for partitioning of the state space into “solve-for” 
and “consider” parameters, allowed for differences between the formal values and the true val- 
ues of the measurement noise, process noise, and a priori solve-for and consider covariances, 
and explicitly partitioned the errors into subspaces containing only the influence of the mea- 
surement noise, process noise, and a priori solve-for and consider covariances. In this work, 
we explicitly add sensitivity analysis to this prior work, and relax an implicit assumption that 
the batch estimator’s anchor time occurs prior to the definitive span. We also apply the method 
to an integrated orbit and attitude problem, in which gyro and accelerometer errors, though not 
estimated, influence the orbit determination performance. We illustrate our results using two 
graphical presentations, which we call the "variance sandpile” and the “sensitivity mosaic,” 
and we compare the linear covariance results to confidence intervals associated with ensemble 
statistics from a Monte Carlo analysis. 


NOTATION CONVENTION 


Style 

Example 

Connotation 

Plain 

a, A 

Scalar real number 

Bold 

x, X 

Physical vector, i.e. an “arrow” in 3 -dimensional, physical space 

Calligraphic 

1 ,C 

Coordinate frame in which physical vectors may be expressed 

Sans Serif 

P.0 

Column, row, or matrix of scalars 


INTRODUCTION 

In Refs. 1 and 2, Markley and others proposed a novel approach to linear covariance analysis, for both batch 
and sequential estimators. Unlike the methods described in Maybeck 3 and Gelb, 4 Markley, et al.’s method 
allowed for the partitioning of the state space into “solve-for” and “consider” parameters. Unlike the method 
described in Tapley, et ah, 5 Markley, et al.’s method allowed for differences between the formal values and 
the true values of the measurement noise, process noise, and a priori solve-for and consider covariances. 
Unlike the method described in Jazwinski, 6 the method of Markley et al. explicitly partitions the errors into 
subspaces containing only the influence of the measurement noise, process noise, and a priori solve-for and 
consider covariances. 

In this paper, we review and extend in two directions the results of Refs. 1 and 2, by explicitly adding sen- 
sitivity analysis, and by generalizing the earlier results, which implicitly assumed that the batch estimator’s 
anchor time occurs prior to the definitive span. We also apply the method to an integrated orbit and attitude 
problem, in which gyro and accelerometer errors, though not estimated, influence the orbit determination per- 
formance. We illustrate our results using two graphical presentations, which we call the “variance sandpile” 
and the “sensitivity mosaic,” and we compare the linear covariance results to confidence intervals associated 
with ensemble statistics from a Monte Carlo analysis. 


* Flight Dynamics Analysis Branch, NASA Goddard Space Flight Center, Code 595, Greenbelt, MD 20771. 
+GN&C Systems Branch, NASA Goddard Space Flight Center, Code 591, Greenbelt, MD 20771 
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REVIEW AND UPDATE OF METHOD 


We summarize the generalized linear covariance analysis algorithm of Refs. 1 and 2 here. The algorithm 
performs a general covariance analysis linearized about a given reference trajectory, X(£). It assumes that 
true values and formal values may differ for the a priori covariance, measurement noise covariance, and pro- 
cess noise power spectral density. We denote these parameters as P 0 _, P 0 _, R, R, Q, and Q, respectively. 
The (error-free) dynamics and measurement partials, A (t) and H(f), and the process and measurement noise 
covariances, may be time-varying*. The algorithm assumes that linear, possibly time-varying, transforma- 
tions partition the state space into a “solve-for” state, s (t), which is to be estimated from the measurements, 
and a “consider” state, c(f), which will not be estimated. The mapping of the state-space into solve-for and 
consider subspaces is defined according to 

s (£) = S (f)x(f), c(f) = C(f)x(f) (1) 

where x(f) = X(f) — X(f) is the linear deviation of the estimate from the true state. Mapping back into the 
original state space may be accomplished as follows: 

W) = [ eft] ] ’ M "W = [S(t), C(i)] (2) 

x(t) = S(f)s(f) + C(f)c (t) (3) 

The main computational result of the original algorithm is a partitioning of the contributions to the total error 
from a priori uncertainty, measurement noise, and process noise. We denote the true and formal covariances 
of these error partitions as P 0 , P a , P v , P, , P„., and P w , respectively. The total true variance and total formal 
variance are, respectively, the sums of these partitions: 

P = Pa + Pt) + Pt«) P = P„ + P„ + P„, (4) 

Another way of looking at the results is via difference between true and formal values of the partitions, 
defined as 

AP Q = SP a S T -Pa, AP„ = SP„S T — P„, APu, = SPa,S T — Pu, (5) 

It is perhaps easiest to understand the algorithm as it is applied to the sequential Kalman filter, so we begin 
with this. We then summarize the application of the method to the batch estimator, and generalize the batch 
application so as to accommodate an arbitrary location of the anchor time. We also generalize the batch form 
of the algorithm so that it will propagate the true and estimated covariances, and the sensitivity matrix, over a 
specified time span that does not necessarily correspond to the definitive span. Next, we add to the previous 
work a computation of the sensitivity to a priori errors, and show how the previous work can be used to study 
sensitivity to measurement and process noise. 

Sequential Kalman Filter Form Let 

^ss(ti,t 0 ) = S(fi)<t>(fj,f 0 )S(f 0 ), H s (fj) = H(fj)S(fj) (6) 

and 


ti— l) — 

f <t>ss(ti,T) Q(r)0g S (L,r)dr 

Jti_ 1 

( 7 ) 

ti— l) = 

f 4 , (L,r)Q(r)0 T (£ i ,r)dr 
Jti- 1 

( 8 ) 


*The assumption that A (£) and H (t) are error-free is not as limiting as may at first seem, since use of incorrect parameters or model 
truncations by the estimator will be accommodated in the consider state analysis. 
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where the state transition matrix , t Q ) is derived from Aft) in the usual way, e.g. as specified in Ref. 7. 


Then, propagate the various covariance partitions between measurement updates as follows: 

Pa(tr) = P a (t)") = P 0 _ (9) 

Pa(*D = HUA-JPaittJ&fati-i), P a (fJ") = P 0 _ (10) 

p v (t~) = i)p„(*^Li)<t > ^ JS (*», i), p„(tr) = o (ii) 

p v{ti) = <u(ti,ti-i)p„(t^i)<i> T (ti,ti_i), p,(tr) = o (12) 

Pt«(ij r ) = + Qd(ti,ti-l), p w{ti ) = 0 (13) 

p«x*n = ®{u,u-i)p + Q d (tj,tj_i), p w (tr) = o (i4) 

Note that the a priori covariance is applied only to the a priori partitions, and that process noise is applied 
only to the process noise partitions. 

At each measurement, update the various partitions with a Kalman gain based on the formal covariances: 

Kj = P(i j” ) (i») [H s {U)P(t ~ ) Hg (f j) + R(t,)] _1 (15) 

P a {tt) = [I - KjH s (fj)]P a (f)“)[l - K,H s (f,;)] T (16) 

P a(t+) = [i - S(fj)K i H(f i )]P a (f“)[l - S(t i )K i H(t i )] T (17) 

P„(f+) = [i - K i H s (fj)]P t) (f)“)[l - K i H s (t i )] T + KiR(ti)KT (18) 

P v(t+) = [I - S(f,)KjH(f i )]P t ,(f“)[l - S(fj)K,;H(£j)] T + S(£ i )K l R(£ i )KjS T (£ i ) (19) 

P w (4) = [I - KiH a {ti)]Pw{ti)[\ - KjH s (£j)] T (20) 

Pw(4) = [I - S(t i )K i H(^)]P-(^)[l - S(t i )K i H(t i )] T (2D 


Note that the measurement noise is applied only to the measurement noise partitions. 

Batch Estimator Form The batch update maps all the innovations to an anchor time, which is associated 
with any a priori information, but which is not necessarily prior to the definitive data span; to emphasize 
this, we denote the anchor time with the subscript V rather than ‘o’ The a priori and measurement noise 
partitions can be expressed as follows: 


K, = [P.-I + £ lU/jch.j/,./,.); ^ ss ( tl ,umu)R- 


( 22 ) 


p a{ti) = [I - E i K<H i (t i )<l>„(t i ,t,)]P._[l - EiKiH a (ti)<t> ss (ti,U)] T (23) 

Pa(^) = [I - E i S(ti)K i H(t i )<l>(t i ,t.)]P.-[l - (24) 

Pv(t+) = EiKiR(t<)KT, P„(i+) = ^S^R^KJS^) (25) 


The contribution to the batch estimator’s formal variance from process noise is zero, since the batch filter 
does not model process noise; however, if there is really process noise present, then the true variance due to 
process noise is 


p w (ti) = t 


t 2 , ti) Qd{t*;t 2 ,t 2 ) 



(26) 


where 

and 


T = 


S^OKxH (h)Hh,U) S(£ 2 )K 2 H(£ 2 )ch(£ 2 ,£*) 


Q(t (£* i i'i , tj ) 


jmniiuij) t)Q(t )<S> T ( tj, T )dr f* < U,U < tj, 
/max(t j)tj ) c t , (ii,T)Q(r)ct> T (£ i ,r)dr U < U,tj < f*, 

0 otherwise 


(27) 

(28) 
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which reduces to 


Qd(tj,U)<t> T (tj,tj) U <U< tj, 

<K ti,tj ) Q d (tj , f * ) f * < tj < U , 

Q d(t*;U,tj) = < <t>{ti,tj)Q d (t*,tj) U < tj < f*, (29) 

tj < U < f*, 

0 otherwise 

where we note that if tj - tj > £*, then the definition above reduces to the usual definition of the process 
noise covariance, i.e. Q d (f*; tj, tj) = Qd(tj, £*). If t, = tj < f* then Q d (f*; tj, tj) = Q d(t*,tj), i.e. we have 
the usual process noise covariance, but it increases as time flows backwards. In the same fashion, the first and 
fourth cases in Eq. 29 produce the same numerical result, with the former corresponding to forward and the 
latter to reverse time flow. The same holds for the second and third cases in Eq. 29. 

To propagate the estimated and tme covariances, respectively, over a specified time interval (forward or 
backward in time), use the following: 

P a(t) = <t> ss (t,t*)P a (t+)Q T ss (t,t*), P a (t) = <t>(t,t*)P a (t+)<t> T (t,t*) (30) 

P„(f) = P v(t) = Ht,t*)P v {ti)<S> T (t,U) (31) 

P w {t) = 4>(t,U)P w (tt)4> T (t,U) + + N T d (t)0 T (t,U) + Q d (U;t,t) (32) 

where 

N d (t) = -52jS(tj)KjH(tj)Qd(t*;t,tj) (33) 

We note that when we allow for time reversal here, we are considering the postdiction problem; the fact that 
the process noise terms increase the uncertainty as time runs backwards reflects our increasing uncertainty 
as we propagate further into the past. One should not infer from this work that we take any position as to 
whether or not a physical Brownian process is diffusing backward in time. 

Sensitivity Matrix Time and Measurement Update The sensitivity matrix shows the linear sensitivity of 
the solution at a specified time to mismodeling of the distributions of the a priori parameters (both solve-fors 
and considers). For the sequential filter, it is given by 

£ a {U ) = [I -S(£ i )K i H(£ i )]<l>(£ i ,£ i _ 1 )Z 0 (£ i _ 1 ), E 0 (f 0 ) = [S(f 0 ), C(f c )] (34) 

and for the batch estimator, it is 

5=„(t) = S(t)0(t, t*)[l- Et 1 S(*i)KiH(t i )d>(^ ) £*)][S(t*),C(£*)] (35) 

Although it is possible to compute the sensitivities to each particular measurement and process noise sam- 
plc' , this does not appear to be particularly useful. Instead, suppose R and R have the same structure, but 
differ only by a scalar multiplier, i.e., R = r ■ R and R = r ■ R (for example, let R = I). Then, for the sequential 
filter, 

AP„ = SP t ,S T - P„ (36) 

will contain terms of the form KK T (r — f). In this case, if one chooses (r — f) = 1, then AP„ will represent 
the sensitivity to measurement noise mismodeling. Similarly, when Q = q ■ I, and Q = q ■ I, then (q — q) 
will factor out of the process noise partition, so that if one chooses q — q — 1, then AP„ will represent the 
sensitivity to process noise mismodeling. For the batch estimator, 

AP„ = KRK T (r-f). (37) 

In this case, if one chooses (r — f) = 1, AP,, will represent the sensitivity to measurement noise mismodeling 

across the entire definitive data span. Similarly, when Q = q ■ Q, q will factor out of the (true) process noise 
partition, so that if one chooses q = 1, AP^, (which equals P„, since the batch has no process noise) will 
represent the sensitivity to process noise mismodeling across the entire interval of interest. 

I It is particularly easy to show that for the batch estimator, the sensitivity to measurement and process noise inputs is 
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INERTIAL NAVIGATION SYSTEM MODEL 


Next, we apply the generalized method to a study of the effect of inertial navigation system (INS) parame- 
ters on an entry trajectory. This section describes the model of the INS. As part of this description, we address 
some issues related to attitude error modeling which appear to us to have never been adequately described in 
other works. 

Our model includes additive bias and noise terms, and multiplicative errors. One of the multiplicative 
errors is associated with uncertainties in the scale factors used in the devices to convert their internal units 
into output units. The other error considered is that caused by small misalignments in the mounting of the 
gyros and accelerometers, so that they are not quite orthogonal to one another. This model shares many of 
the characteristics of the model used in Ref. 8, where Space Shuttle navigation performance with strapdown 
systems was evaluated, as well as with the model described in Ref. 9, a popular short course on inertial 
navigation systems. Similar models for platform-type inertial systems are described in Refs. 10 and 11. 
A more complicated model that shares some similarities is that of Ref. 12, in which a Kalman filter for 
calibration and alignment of inertial navigation systems is described. 

Gyro Model 

INS are generally of two types, “platform” and “strapdown,” which use the gyros in a somewhat different 
fashion. The platform type consists of a set of gimbals that interconnect the vehicle body to a rigid platform 
on which the gyros are mounted. This arrangement allows the axes along which the gyros are sensitive to 
rotate freely with respect to the vehicle body. A feedback system attempts to null the gyro outputs, which 
thus maintains the gyro platform approximately fixed in inertial space. The vehicle attitude with respect to 
inertial space is then given by reading out the gimbal angles. The strapdown type, which has largely replaced 
the platform type in current usage, dispenses with the gimbal system and rigidly fixes the gyros to the vehicle 
body. The strapdown INS computationally integrates the gyro outputs so as to keep track of the orientation 
of the body with respect to the inertial frame. 

In each type of INS, errors in the gyro output affect the determination of the vehicle attitude, but in a 
somewhat different manner. In a platform INS, gyro output errors cause the platform to physically drift with 
respect to the true inertial frame. The strapdown INS measures at each instant an incremental change in the 
angular position of the body with respect to some inertial frame. Thus, errors in the gyro outputs create error 
in the computational accumulation of these angular velocity increments. 

To illustrate the implications of this difference, we assume (without loss of generality) that the INS rep- 
resents the attitude as a direction cosine matrix that maps from the INS case framed to the inertial frame, 
and we denote its estimate of the attitude as M x . One may view this estimate as composed of either of the 
following successive rotations: 

m|=m£m£, or m| = MfM J (38) 

We assume that we can represent and Mf as small angle rotations: 

M? = I — 9 X and M| = I - ip x (39) 

where the superscript ‘ x ’ indicates that the elements of the vector are arranged in a skew-symmetric matrix 
such that a x b = a x b. Note that 9 and xjj are in general two different sets of small angle rotations - it is 
not necessarily the case that they are the same vector expressed in two different coordinate systems^. The 
(coordinate-independent) derivatives of these vectors with respect to the case and inertial frames are related 

by i c i c 

— 9 = — 0 + x w c x 9 and — r/> = — xp + I u: c x (40) 

dt dt di di 

•I- For either type of INS, the case frame is rigidly fixed to the vehicle body, but may not be aligned with the body axes. 

§ However, if these rotations arise from two optimal attitude estimators processing the same data from the same initial conditions, the 
constraint that both estimators produce the optimal estimate implies that 0 and ip do become the same vector expressed in two different 
coordinate systems. 
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where we denote the angular velocity of frame C in frame 1 as x uf . Similarly, the angular velocity estimate 
of the INS, x u3 c , will be corrupted by the gyro output error, and we may view this estimate as either 

x u c = x uj c + c lo c or x uj c = x u x + x af (41) 

Thus, the coordinate-independent angular velocity error vector, S X U) C = x ui c — x u> c , may be viewed as either 

S x u: c = c uj c or 8 x oj c = x lo x (42) 


Recalling the description of how the platform INS uses the gyro output, it is clear that for the platform 
system, the gyro error corresponds to a drift of the platform with respect to inertial, and since the platform 
represents the INS’s knowledge of the inertial frame, it is clear that the platform’s drift corresponds to x lj x . 
Thus, ip must be related to x aj x , according to 


= V (43) 

at 

Recalling that the strapdown INS integrates the case-fixed gyro outputs, it is clear that a constant gyro drift 
from a strapdown gyro will appear to be constant to an observer who is also fixed in the case frame, and hence 
gyro drift from a strapdown INS corresponds to c ui c , or 

C d n n 

-rd = C ^ C (44) 

at 

Thus, depending on whether the INS is of the platform type or the strapdown type, a constant gyro output bias 
will appear to be constant relative to the inertial frame, or constant relative to the case frame, respectively. 
These gyro biases will integrate into a misalignment ip with respect to the inertial frame (platform INS), or 
into a possibly different misalignment 6 with respect to the case frame (strapdown INS). Note that one is free 
to write the equations of motion with respect to any frame one chooses, so that for example the following are 
completely permissible expressions: 

— 0 = c af + x l o c x 9 and —ip = x u x — x uj c x ip (45) 

at at 

Henceforth, we only consider strapdown INS. For many such systems, the gyro error may be modeled as 

PJt ~ C 4 — t*gC + Sg X u>£ + r x aJ^ + , (46) 

where b g c is an angular velocity in the case frame that biases the gyro, S, ; is the gyro scale factor matrix, Y g 
is the gyro nonorthogonality matrix, x ui£ is the representation of x u £ in coordinates fixed to the INS case, 
is a Gaussian white noise that produces angular random walk, and where 



Sgx 

0 

0 


0 

0 

o ' 

S s = 

0 

Sgy 

0 

and f g = 

~ 'Ygyz 

0 

0 


0 

0 

N 


. Igzy 

/ " Ygzx 

0 


Here, the elements of the matrix S g represent scale factor errors, and the elements of the matrix T, ; represent 
nonorthogonality errors. The latter arise from the fact that (for example) the y gyro will measure not only the 
angular velocity about C y as it is intended to, but also a small projection of the angular velocity about C x \ 

€ Here we assume that as part of the INS calibration procedure, the x gyro is defined to be co-linear with the case x axis, the x-y 
case frame is defined to include the y gyro, and misalignments of the y and z gyros from the two axes orthogonal to this reference are 
measured. In practice, the case frame is defined by an alignment device such as an optical comer cube attached to the outside of the 
case, and the relationship between the frame defined externally by the optical cube and the internally-defined case frame is loaded into 
the INS's firmware. 
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Accelerometer Model 

Since we have assumed a strapdown INS, the accelerometers measure the case coordinates of the “sensed 
acceleration”, a c, i.e. that part of the second time derivative of r, the INS position vector, with respect to an 
inertial frame, which is due to non-gravitational, or “contact,” forces. If the INS is located at the spacecraft 
center of gravity (c.g.), r represents the spacecraft position vector as well. If not, the accelerometer will 
measure not only the external acceleration on the spacecraft, but also the linear acceleration at the c.g. to 
INS moment arm due to any spacecraft rotational motion. This must be subtracted from the INS measured 
acceleration, typically by the navigation software. This model does not explicitly include the effect INS errors 
have on this c.g. offset calculation. 

For many strapdown accelerometers, we may assume that the error may be modeled as 


8&c — t *aC + S a ae + V a a.Q 1 


(48) 


where b a c is the accelerometer bias, S 0 is the accelerometer scale factor matrix, f a is the accelerometer 
nonorthogonality matrix, and where 



5ai 

0 

0 


0 

0 

0 ' 

S a = 

0 

S ay 

0 

, r a = 

'Yayz 

0 

0 


0 

0 

S az 


' lazy 

lazx 

0 


Combined Model 

It is convenient to recast Eqs. 46 and 48 as follows. Define D( x u;g ) = diag( x CL>g), and s g = [ S gx , S gy , S^] 7 . 
Then S 9 x u;g = D( x wg)s g . Similarly, define D(a c ) = diag(a c ), and s a = [S ax , S ay , Sazf- Then S a a c = 
D(a c )s a . Define j g and F( x u>g) such that T 9 x u>g = F( x u>g) 7 9 , i.e. = [j gyz , ^ gzy , 'y gzx ] T , and 


F( T «g) 


0 0 0 

-ML o o 

0 x u4 - x u,g ; 


(50) 


Similarly, define and F(ac) such that T 0 ac = F(ac)7 a - Then, Eqs. 46 and 48 may be written as 

8 x u£ = b gC + D( x u>g)s g + F( x wg) 7 9 + (51) 

<5a c = b oC + D(a c )s a + F(a c ) 7 a (52) 


Now, let x fl = [9 J C , b^ c , sj, 7 9 ] T - Then, one can rewrite Eq. 46 as 

( C i«c \ 


di b 9C 


\ J 


0 3x3 Is D( x u;g) F( x u;g) 

O 9 X 12 


I 3 

09x3 


*9 — Ag x g + B g s u , 


(53) 

(54) 


where we have chosen to write the equations of motion for the gyro misalignment and drift with respect to 
the case frame, and we recall that since the scale factor and nonorthogonality coefficients are only shown as 
column arrays for convenience, they are not properly viewed as physical vectors that take derivatives in any 
particular coordinate system. In similar fashion, let x a = [b^ c , s^, 7 ^] T , with x a = ^x a = 0, and cast the 
acceleration error <5a c as 


<5a c = [ l 3 D(a c ) F(a c ) ] x a 
= H a x a . 


(55) 

(56) 
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We refer to the combination of the two models above as our Inertial Measurement Unit (IMU) model. 
When we integrate these results to produce position and velocity errors, we have a full INS model. When 
an INS incorporates sensed acceleration into its computation of the total acceleration on the vehicle, its gyro 
misalignment will corrupt the calculation: 

-jtVj = - ^fx + M x c k c (57) 

at r 6 

« -^r x + G{r I )5r I + Mj[l - ](a c + <5a c ) (58) 

where the gravity gradient matrix, G(rx), is the Jacobian of the two-body gravity term, —firj/r 3 . Thus, to 
first order, we may approximate the total acceleration error equations of motion as 


T d 

— dvx = G(rx)c>rx + Mja ^G c + M ^5a. c 
These models may be combined as follows: 


/ 



03x3 

I 3 

03x9 

03x12 


( \ 




S^x 

x a 


G(rx) 

03x3 

M?H a 

MjXag 


Svz 

_|_ 

Ol5x3 



Ol2x3 

Ol2x3 

09x9 

Ol2xl2 


Xa 

1 

B 9 . 

V 

*, J 


Oi2x3 

Oi2x3 

09x9 

A 3 . 


V x s / 




x = Ax + Be w 


where X ag = [a£ , O 3 X 9 ] couples the gyro and accelerometer errors. 


(59) 


(60) 

(61) 


EXAMPLE 

Next, we describe our example problem, then use the generalized covariance method to analyze the per- 
formance and sensitivities of the inertial navigation system. 


Description of Example Problem 

In our example, we assume that a global navigation satellite system type of receiver provides vector mea- 
surements of position throughout an entry trajectory, with measurement noise standard deviation three times 
less than the estimator assumes. Figs. 1 and 2 show some salient characteristics of the example. The sim- 
ulation begins after the de-orbit burn, and one sees in the top plot of Fig. 1 that the vehicle reaches a peak 
deceleration of about 12 x Earth’s surface gravity about half-way through the entry. The vehicle executes 14 
bank maneuvers, which are a combination of rolls and yaws, as the top and bottom plots of Fig. 2 show, so 
as to control its landing point. The middle and bottom plots of Fig. 1 show accelerations from side forces 
created by the vehicle’s lift during these bank maneuvers. 

A Monte Carlo simulation produced the gyro and accelerometer error data that Figs. 1 and 2 also show. 
This simulation assumed that the elements of 9c(t 0 ), b g c, s g , 7 g , b a c, s a , ~f a and e UJ are zero mean, timewise 
uncorrelated, Gaussian-distributed random processes, and are uncorrelated with one another. Table 1 lists the 
relevant simulation parameters that we used. 


Navigation Analysis of Sample Problem 

In using any linear covariance analysis method, we always recommend at least one actual simulation as a 
check against linearization problems. In the present work, we generate random deviations from the reference 
as initial conditions for each of 25 Monte Carlo cases. We integrate each deviated case, and use this as truth 
for measurement simulation and estimation error generation. From these, we generate the time series of 
estimation errors and residuals for each case. 



Table 1 Simulation Parameters 


Parameter 

Standard Deviation 

Units 

Gravitational Constant 

4.305 x 10 4 

km 3 /sec 2 

Measurement Time Interval 

2 

sec 

True Position Measurement Noise 

0.305 

m 

Formal Position Measurement Noise 

0.914 

m 

Initial Position Error 

30.5 

meters 

Initial Velocity Error 

3.05 

cm/sec 

Accelerometer Bias 

60 

pg 

Accelerometer Scale Factor 

500 

ppm 

Accelerometer Nonorthogonality 

10 

ppm 

Initial Gyro Angular Error 

42 

arcsec 

Gyro Bias Drift 

0.01 

deg/hr 

Gyro Scale Factor 

33 

ppm 

Gyro Nonorthogonality 

20 

ppm 

Gyro Random Walk Intensity 

0.025 

deg/hr 1 / 2 


Description of Plots We will be discussing a number of plots below, and it will be helpful to describe our 
plot conventions in some detail before delving into the results themselves. 

A typical output of Monte Carlo simulations is the ensemble of the time histories of the true estimation 
errors, along with their formal errors, which are the standard deviations from the diagonals of the estimator’s 
covariance matrix. Fig. 4 is an example of such a result. In this and similar figures to follow, we plot the true 
error time series as a sheaf of yellow lines, and the ensemble mean of the true errors as a single blue line. 
We plot the empirical ensemble deviation of the true errors as a thick cyan band that covers ±1 to ±3cr’s. 
We show twice the true standard deviation from the linear covariance analysis as a series of black + marks. 
We show twice the ensemble mean of the formal standard deviations from the Monte Carlo analysis as a 
series of green + marks. We compute a 95% confidence interval based on a hypothesis test of the variance 
of the Monte Carlo results against the true variance from the linear covariance analysis, which we show as 
a translucent red band bounded by red + marks. These sorts of results are available from existing linear 
covariance and Monte Carlo simulation procedures. 

We have chosen to illustrate the capability of the generalized covariance method to partition the covari- 
ance using what we call “variance sandpiles,” which are stacked area charts showing the time series of each 
solve-for variance’s contribution from a priori error variance, measurement noise variance, and process noise 
variance. Fig. 7 is an example of such a result. Here, we plot the components of the true variance as a posi- 
tive sandpile, and the components of the formal variance as an inverted sandpile on the negative ordinate, and 
relabel the negative ordinate to indicate this I . 

A useful result of many types of covariance analysis is an error budget showing the sensitivity of each 
solve-for parameter to errors in the consider parameters. We show such results as “sensitivity mosaics,” 
which are checkerboard plots of the log magnitude of the elements of the sensitivity matrices. Fig. 8 is an 
example of such a plot. 

Results We performed a number of studies of our example problem using the generalized method, which 
we briefly summarize before discussing our final results. First, we solved for all of the states in the INS model 
to study the observability of the various states. From this work, we concluded that accelerometer bias was 
weakly observable and accelerometer scale factor was strongly observable along the two directions generally 

I Depending on the relative magnitudes of these results, we could also illustrate the sandpiles in other ways, e.g. if all the delta 
variances were positive, the sandpile could show the formal variance, and the deltas due to each component. If all the deltas were 
negative, the sandpile could show the true variance, and negatives of the deltas. 
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normal to the drag acceleration. Accelerometer nonorthogonality and all of the gyro rate errors were not 
observable. Gyro misalignment was strongly observable, along with position and velocity. Based on these 
results, we next examined various combinations of observable states as solve-for parameters. Among these 
combinations, we studied (1) solving for position, velocity, accelerometer bias and scale factor, and gyro 
misalignment; (2) solving for only position and velocity; and (3) solving for linear combinations of the 
observable accelerometer and gyro states. The latter was particularly easy with the generalized formulation, 
since we needed only to modify the S and C matrices. From this work, we concluded that the best performance 
was available from case (1). Figs. 3-8 show the results from this case we obtained with the sequential form 
of the method. Finally, we studied this case using the batch estimator. Although this produced comparable 
results, we found that the large size of the T matrix, by which the batch form captures the effect of process 
noise, caused the batch analysis to take considerably longer to run. 

Discussion There are several notable features of Fig. 3, which shows the Monte Carlo and linear covari- 
ance results for the position states**. The agreement between the true standard deviation from the linear co- 
variance analysis (black + marks) and the ensemble mean of the formal standard deviations from the Monte 
Carlo analysis (green + marks) indicates that nonlinearities are not significantly affecting the results. These 
are bounded by the Monte Carlo confidence interval (red band), indicating that the null hypothesis, which is 
that the Monte Carlo cases actually are governed by the moments of the linear covariance analysis, cannot 
be rejected at the 95% confidence level. The large sawtooth increases in the errors between measurement 
updates lead us to suspect that process noise, which only comes from gyro random walk in this example, is a 
strong factor influencing performance. 

The position variance sandpiles (left column of Fig. 6) allow us to confirm that the spikes are due to 
process noise. They also show that measurement noise is a significant contributor to the formal variance, but 
not the true variance (recall that the simulation was run with three times less measurement noise on the true 
measurements than the filter assumed). Referring to the sensitivity mosaic in Fig. 8, we can see that the IMU 
parameters, especially gyro drift, are significant contributors to position and velocity errors. 

Considering now the accelerometer and gyro errors in Figs. 4 - 5, we see that the ^-component of the 
accelerometer scale factor error somewhat exceeds the 95% confidence limit during the initial convergence 
of the filter, possibly indicating stronger nonlinear effects not captured by the linear covariance analysis 
(Fig. 4). Once the ^-component of the gyro angular error converges after t = 40 (Fig. 5), the confidence 
interval violation for the accelerometer scale vector stops, suggesting that the cross-product between gyro 
angular error and sensed acceleration is playing a role. The sandpiles for the weakly observable accelerometer 
parameter variances (Fig. 7) are dominated by a priori error, while the strongly observable gyro angular error 
variance sandpile (right column of Fig. 6) is dominated by process and measurement noise, as with the other 
strongly observable states. 


SUMMARY AND CONCLUSION 

We have updated and extended previous work on linear covariance analysis and used our generalized 
method to study the navigation performance of an inertial navigation system dining an entry trajectory. As 
part of this work, we clarified some issues concerning the equations of motion for modeling attitude errors. 

We showed that sensing drag acceleration during the entry, along with position measurements from a GPS- 
like device, allows the INS to update its attitude error states. We found that the linear covariance analysis 
results agreed well with a 25-case Monte Carlo simulation, indicating the lack of any strongly nonlinear 
effects in our example problem. Using the unique capability of the generalized method to partition the errors, 
we were able to quantify the extent to which process and measurement noise dominate the performance in 
estimating the strongly observable states, and to quantify the extent to which a priori error strongly affects the 
performance in estimating the weakly observable states. We illustrated the sensitivity of the INS to various 
unobservable parameters, and showed that position and velocity errors are strongly sensitive to gyro bias. 


** Velocity results are generally similar to position, so we omit them for brevity. 
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Figure 1 Monte carlo simulation of case coordinates of accelerometer errors (left scale, solid 
blue lines) and their formal 3<r values (left scale, dotted blue lines), shown with reference accel- 
erations (right scale, dashed red lines); the x, y, and z coordinates are shown in the top, middle, 
and bottom subplots, respectively. 


12 


Earth G's 


arcsec/sec arcsec/sec arcsec/sec 


IMU Model Gyro Rate Errors 


60 


40 


20 - 


Body Rate Errors [-] (30 cases); Lincov 3o Reference Body Rates [- 


I 


-■ - ' l 

, / ' 


t i r 


> 


T i - 

'I 

l| 







■H — 


— 1 40 
- 20 
- 0 
- -20 

- -40 

- -60 
- -80 
- -100 


-20 


20 


40 


60 


80 


100 


120 


140 


160 180 


-120 




Elapsed Time [sec] 

Figure 2 Monte carlo simulation of case coordinates of gyro rate errors (left scale, solid blue 
lines) and their formal 3<r values (left scale, dotted blue lines), shown with reference body rates 
(right scale, dashed red lines); the x, y, and a coordinates are shown in the top, middle, and 
bottom subplots, respectively. 
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Figure 3 Position Errors from 25-case Monte-Carlo Simulation. 
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Figure 4 Accelerometer Scale Factor Errors from 25-case Monte-Carlo Simulation. 
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Figure 5 Gyro Angular Errors from 25-case Monte-Carlo Simulation. 
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(c) Y -component of Inertial Position Error 
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(e) Z-component of Inertial Position Error 


(f) Z-component of Case-Fixed Gyro Angular Error 


Figure 6 Position Error (a,c,e) and Gyro Angular Error (b,d,f) Variance Sandpiles (plotted vs. sample number). 
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(a) X -component of Case-Fixed Accelerometer Bias 


(b) X-component of Case-Fixed Accelerometer Scale Factor 




(c) Y -component of Case-Fixed Accelerometer Bias 
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Figure 7 Accelerometer Bias (a,c,e) and Scale Factor (b,d»f) Variance Sandpiles (plotted vs. sample number). 
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Figure 8 Sensitivity of solve-for errors to variations in a priori parameter specifications. The 
solve-for indices are as follows: 1-3, position; 4-6, velocity; 7-9, accelerometer bias; 10-12, 
accelerometer scale factor; 13-15, gyro angular error. 
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